# cython: boundscheck=False
# cython: wraparound=False
from cython cimport boundscheck,wraparound

import numpy as np
from libc cimport math
from libc.stdlib cimport malloc


DEF PI = 3.141592653589793
DEF M_3_2_PI=4.71238898038
DEF M_2__PI= 6.28318530718
DEF NOTDEF=-1024
DEF INF=1.7976931348623158e+308
DEF density_thresh=0.7

cdef struct coorlist:
    int x
    int y
    void *next

cdef class Rect:
    cdef public double x1,y1           #矩形主轴反方向上的端点
    cdef public double x2,y2           #矩形主轴正方向上的端点
    cdef public double cx,cy           #矩形的中心点
    cdef public double width               #宽度
    cdef public double theta        #矩形的方向（弧度）
    cdef public double precision
    cdef public double p
    cdef public double dx,dy



"""
double[:]指示kernel为typed memoryview
"""
def gaussian_kernel(int size,double[:] kernel,double sigma,double mean):
    cdef double sum = 0.0,sum_div
    cdef unsigned int i
    for i in range(size):
        val = (i - mean)/sigma
        kernel[i] = math.exp(-0.5 * val ** 2)
        sum = sum + kernel[i]
    if (sum >= 0.0):
        for i in range(size):
            kernel[i] = kernel[i]/sum
    return


def gaussian_sampler(double[:,:] image, double scale,double sigma_scale):
    cdef:
        double prec,mean,xx,yy,sum
        int N,M,ysize,xsize, x,y,i,j,xc,yc
        int h,n,double_x_size,double_y_size
        double[:] kernel
        double[:,:] aux
        double[:,:] out
        #double *aux_arr
        #double *out_arr
    ysize,xsize=image.shape[:2]
    N=int(math.ceil(xsize*scale))
    M=int(math.ceil(ysize*scale))
    aux=np.zeros([ysize,N],dtype=np.double)
    #aux_arr=<double*>malloc(ysize*N*sizeof(double))
    #aux=<double[:ysize, :N]>aux_arr
    out=np.zeros([M,N],dtype=np.double)
    #out_arr=<double*>malloc(M*N*sizeof(double))
    #out=<double[:M,:N]>out_arr
    if scale<1.0:
        sigma=sigma_scale/scale
    else:
        sigma=sigma_scale
    prec=3.0
    h=int(math.ceil(sigma*math.sqrt(2.0*prec*math.log(10)))  )
    n=1+2*int(h);    #高斯核的尺寸
    kernel = np.zeros(n,dtype=np.double)

    #####先对X轴减采样
    double_x_size=2*xsize
    double_y_size=2*ysize
    for x in range(N):      #x为新图像中的x坐标
        xx=x/scale  #与x对应的原图像中的坐标
        xc=int(math.floor(xx+0.5))     #对xx的四舍五入,核操作的中心像素
        ##生成高斯核
        mean= h+xx-xc;
        gaussian_kernel(n,kernel,sigma,mean)
        for y in range(ysize):
            sum=0.0
            for i in range(n):  #以xc为中心的像素与核进行运算
                j=xc-h+i;
                while j<0: j=j+double_x_size
                while j>=double_x_size: j=double_x_size-j
                if j>=xsize: j=double_x_size-1-j
                sum=sum+image[y,j]*kernel[i]
            aux[y,x]=sum
    ###经过aux之后得到ysize*N的图像，X轴方向减采样了
    ######对y轴减采样
    for y in range(M):
        yy=y/scale
        yc=int(math.floor(yy+0.5))
        ##生成高斯核
        mean = h + yy - yc;
        gaussian_kernel(n,kernel,sigma, mean)
        for x in range(N):
            sum=0.0
            for i in range(n):
                j=yc-h+i
                while j<0: j=j+double_y_size
                while j>=double_y_size: j=double_y_size-j
                if j>=ysize: j=double_y_size-1-j
                sum=sum+aux[j,x]*kernel[i]
            out[y,x]=sum
    return np.asarray(out)

def levelline_angle(double[:,:] image,double threshold,int n_bins):
    """
     计算图像的梯度，等高线方向，以及进行梯度伪序操作
     threshold及论文中的rho
     返回值：梯度幅值，梯度方向，梯度伪序排序的像素坐标列表
    """
    cdef:
        int r,c,height,width,i,list_count
        double com1,com2,gr,gc,max_grad,norm
        double [:,:] mag
        double [:,:] ll_angles
        int [:,:] list_p
        coorlist **range_l_s
        coorlist **range_l_e
        coorlist *start
        coorlist *end
        coorlist *list

    ##计算梯度幅值，使用2x2窗口
    print("compute grad.")
    height,width=image.shape[:2]
    #out_image=np.zeros((height,width))  #输出图像
    mag=np.zeros((height,width),dtype=np.double)    #梯度幅值
    ll_angles=np.ones((height,width),dtype=np.double)*NOTDEF   #等高线方向的角度（弧度）
    max_grad=0.0
    ##计算梯度，角度
    for r in range(height-1):
        for c in range(width-1):
            com1=image[r+1,c+1]-image[r,c]
            com2=image[r,c+1]-image[r+1,c]
            gr=com1-com2;   gc=com1+com2
            #计算梯度幅值
            norm=gr*gr+gc*gc
            norm=math.sqrt(norm/4)
            mag[r,c]=norm
            #计算梯度角度
            if norm>threshold:
                ll_angles[r,c]=math.atan2(gc,-gr)
                if norm>max_grad: max_grad=norm
            else:
                ll_angles[r,c]=NOTDEF    #小于阀值的像素角度设为未定义

    ##得到梯度伪序
    #先创建列表
    print("fill bins.")
    list=<coorlist *> malloc(height*width*sizeof(coorlist))
    range_l_s=<coorlist **>malloc(n_bins*sizeof(coorlist *))
    range_l_e=<coorlist **>malloc(n_bins*sizeof(coorlist *))
    if list==NULL or range_l_s==NULL or range_l_e==NULL:
        raise MemoryError()
    for i in range(n_bins):
        range_l_s[i]=NULL
        range_l_e[i]=NULL

    #按梯度值将像素分到对应的bin中
    print("fill list_p.")
    list_count=0
    for c in range(width-1):    #图像的最后一行和最后一列去掉
        for r in range(height-1):
            norm=mag[r,c]
            i=<unsigned int>(norm*n_bins/max_grad)     #计算像素梯度值对应的bin
            if i>=n_bins: i=n_bins-1
            #从链表中取一个元素,放到对应的bin中
            if range_l_e[i]==NULL:              ##当前bin为空时
                range_l_s[i]=&list[list_count]
                range_l_e[i]=&list[list_count]
            elif range_l_s[i]==range_l_e[i]:    ##当前bin有一个元素
                range_l_e[i]=&list[list_count]
                range_l_s[i].next=<void*>&list[list_count]
            else:                               ##当前bin有超过一个元素
                range_l_e[i].next=<void*>&list[list_count]
                range_l_e[i] =&list[list_count]
            list_count=list_count+1
            range_l_e[i].x=c
            range_l_e[i].y=r
            range_l_e[i].next=NULL

    #按bin从高到低，将所有元素排成一列
    list_p=np.zeros([(width-1)*(height-1),2], dtype=np.int)
    list_count=0
    for i in range(n_bins-1,-1,-1):
        start=range_l_s[i]
        while start!=NULL:
            list_p[list_count,0]=start.y
            list_p[list_count,1]=start.x
            list_count=list_count+1
            start=<coorlist*>start.next

    return np.asarray(mag),np.asarray(ll_angles),np.asarray(list_p)



def isaligned(int x,int y,double angle,double[:,:] ll_angles,double prec):
    """
    判断(y,x)像素与给定区域角度值是否一致
    """
    cdef double a
    a=ll_angles[y,x]
    if a== NOTDEF: return False
    angle=angle-a
    if angle<0: angle = -angle
    if angle>M_3_2_PI:
        angle=angle-M_2__PI
        if angle<0: angle=-angle
    return angle<=prec

def region_grow(int x,int y,double[:,:] ll_angles, int[:,:] reg, int[:,:] used, double prec):
    """
    区域生长算法，以(x,y)点开始进行区域生长
    ll_angles:等高线角度图，used:用于标记像素是否使用
    返回值：区域像素列表，区域角度
    同时会修改used
    input parameters：ll_angles is (M,N) level-line angles，
                    reg is a (N*M,2)region point buffer
    """
    cdef:
        int height,width,px_count,i, cx,cy,xx,yy
        double reg_angle, sum_sin,sum_cos
    height,width=ll_angles.shape[:2]
    reg[0,0]=y
    reg[0,1]=x
    px_count=1  #used for region pixel count
    used[y,x]=1
    reg_angle=ll_angles[y,x]  #标量，记录整个区域的角度
    sum_sin=math.sin(reg_angle)
    sum_cos=math.cos(reg_angle)
    i=0
    while  i<px_count:
        #取一个元素进行区域操作
        cx=reg[i,1]; cy=reg[i,0]
        xx=cx-1
        while xx<=cx+1:
            yy=cy-1
            while yy<=cy+1:
                if yy>=0 and yy<height and xx>=0 and xx<width and (not used[yy,xx]) and   \
                    isaligned(xx,yy,reg_angle,ll_angles,prec):
                    ##添加像素点到区域
                    used[yy,xx]=1;
                    reg[px_count,0]=yy
                    reg[px_count,1]=xx
                    px_count=px_count+1
                    ##更新区域的角度
                    sum_sin=sum_sin+math.sin(ll_angles[yy,xx])
                    sum_cos=sum_cos+math.cos(ll_angles[yy,xx])
                    reg_angle=math.atan2(sum_sin,sum_cos)
                yy=yy+1
            xx=xx+1
        i=i+1
    return reg_angle, px_count

def region2rect(Rect rect,int[:,:] region,int px_count,double reg_angle,double precison,double p, double[:,:] grad_mag):
    """
    为区域生成矩形
    reggion:区域像素列表  reg_angle:区域角度
    prec: 精度(弧度)    p: 像素是对齐点的概率，
    返回值：矩形的中心点，方向，宽度，长度
    """
    cdef:
        double sum_x,sum_y,sum,Ixx,nIxy,Iyy,w_sum,weight
        double cx,cy,dx,dy,theta,lambd,angle_diff,l_min,w_min,l_max,w_max,l,w
        int i,x,y

    ##计算区域质心，作为矩形中心
    sum_x=0.0; sum_y=0.0; sum=0.0
    for i in  range(px_count):
        x=region[i,1]; y=region[i,0]
        weight=grad_mag[y,x]
        sum_x=sum_x+x*weight
        sum_y=sum_y+y*weight
        sum=sum+weight
    cx=sum_x/sum
    cy=sum_y/sum
    ##计算矩形方向
    Ixx=0.0; nIxy=0.0; Iyy=0.0; w_sum=0.0
    for i in  range(px_count):
        x=region[i,1]; y=region[i,0]
        weight=grad_mag[y,x]
        w_sum=w_sum+weight
        Ixx=Ixx+(x-cx)*(x-cx)*weight
        Iyy=Iyy+(y-cy)*(y-cy)*weight
        nIxy=nIxy-(x-cx)*(y-cy)*weight       #在惯量矩阵中，反对角线上是-Ixy
    Ixx=Ixx/w_sum; Iyy=Iyy/w_sum; nIxy=nIxy/w_sum
    lambd=0.5*(Ixx+Iyy-math.sqrt((Ixx-Iyy)*(Ixx-Iyy)+4.0*nIxy*nIxy) )
    if abs(Iyy) >abs(Ixx):
        theta=math.atan2((lambd-Iyy),nIxy)
    else:
        theta=math.atan2(nIxy,lambd-Ixx)
    #检查与区域方向的差值，进行纠正
    angle_diff=theta-reg_angle
    while angle_diff<= -PI: angle_diff=angle_diff+2*PI
    while angle_diff> PI: angle_diff=angle_diff-2*PI
    if abs(angle_diff) >precison: theta=theta+PI
    ##计算矩形宽度
    dx=math.cos(theta)
    dy=math.sin(theta)
    l_min=0; l_max=0; w_min=0; w_max=0
    for i in  range(px_count):
        x=region[i,1]; y=region[i,0]
        l= (x-cx)*dx+ (y-cy)*dy     #像素点在矩形主轴上的投影
        w= (x-cx)*(-dy)+(y-cy)*dx   #像素点在幅轴上的投影
        if l>l_max: l_max=l
        if l<l_min: l_min=l
        if w>w_max: w_max=w
        if w<w_min: w_min=w
    rect.x1=cx+l_min*dx     #与主轴正方向相反的端点
    rect.y1=cy+l_min*dy
    rect.x2=cx+l_max*dx     #与主轴正方向相同的端点
    rect.y2=cy+l_max*dy
    rect.cx=cx
    rect.cy=cy
    rect.theta=theta
    rect.dx=dx; rect.dy=dy
    rect.width=w_max-w_min  #注意，这个矩形是非主轴对称的
    rect.precision=precison
    rect.p= p
    return


def reduce_region_radius(Rect rect,int[:,:] region,int px_count,
                         double[:,:] grad_mag,double[:,:] ll_angles,int[:,:] used,
                         double reg_angle,double precision,double p,density_thresh):
    cdef:
        int xs,ys,i,xx,yy
        double rad1,rad2,rad,dist,density
    xs=region[0,1]
    ys=region[0,0]
    ##计算区域半径
    rad1=math.sqrt((rect.x1-xs)**2+(rect.y1-ys)**2)
    rad2=math.sqrt((rect.x2-xs)**2+(rect.y2-ys)**2)
    if rad1>rad2: rad=rad1
    else: rad=rad2
    while True:
        rad=rad*0.75    #减小半径
        ##从区域中剔除在半径之外的像素点,并设置为未使用
        i=0
        while i<px_count:
            xx=region[i,1]; yy=region[i,0]
            if math.sqrt((xx-xs)**2+(yy-ys)**2)>rad:
                used[yy,xx]=0
                #注意，此时只是用列表尾元素替代当前元素，i值不需要增加
                region[i,0]=region[px_count-1,0]
                region[i,1]=region[px_count-1,1]
                px_count=px_count-1
            else:
                i=i+1
        if px_count<2: return NOTDEF,0,False
        region2rect(rect,region,px_count,reg_angle,precision,p,grad_mag)
        dist=math.sqrt((rect.x1-rect.x2)**2+(rect.y1-rect.y2)**2)
        density=px_count/(dist*rect.width)
        if density>=density_thresh:
            break
    return px_count,reg_angle,True

def angle_diff(double a,double b):
    a=a-b
    while a<=-PI: a=a+M_2__PI
    while a>PI: a=a-M_2__PI
    return a

def refine(Rect rect,int[:,:] region,int px_count,double reg_angle,
            double[:,:] grad_mag,double[:,:]ll_angles,int[:,:] used,
           double precision,double p,double density_thresh):
    """
    首先判断矩形中对齐点密度,达到阀值标准则判定为有效线段，否则进行角度容差和区域的调整使得符合标准
    返回值：region,reg_angle,rect,passed
    """
    cdef:
        int xs,ys,n
        double dist,density,angle_s,angle,ang_d,mean_angle
        double sum,s_sum,tau,length
    dist=math.sqrt((rect.x1-rect.x2)**2+(rect.y1-rect.y2)**2)
    density=px_count/(dist*rect.width)
    #如果满足密度标准，不做处理，返回True
    if density >= density_thresh:
        return px_count,reg_angle,True
    #####若不满足标准，首先尝试减少角度容差
    ##计算新的平均角度和容差
    xs=region[0,1]     #区域种子（seed）点的坐标
    ys=region[0,0]
    angle_s=ll_angles[ys,xs]
    sum=0.0; s_sum=0.0; n=0
    for i in range(px_count):
        xx=region[i,1]; yy=region[i,0]
        used[yy,xx]=0   #标记为未使用
        dist=math.sqrt((xx-xs)**2+(yy-ys)**2)
        if dist< rect.width:  ##使用seed点在小范围邻域，半径为矩形宽度
            angle=ll_angles[yy,xx]
            ang_d=angle_diff(angle,angle_s)
            sum=sum+ang_d
            s_sum=s_sum+ang_d**2
            n=n+1
    mean_angle=sum/n
    tau=2.0*math.sqrt((s_sum-2.0*mean_angle*sum)/n+mean_angle**2)     #两倍标准差
    ##重新进行区域生长
    reg_angle,px_count=region_grow(region[0,1],region[0,0],ll_angles,region,used,tau)
    if px_count<2:
        return NOTDEF,0,False
    ##重新计算矩形
    region2rect(rect,region,px_count,reg_angle,precision,p,grad_mag)
    #重新计算密度
    length=np.sqrt((rect.x1-rect.x2)**2+(rect.y1-rect.y2)**2)
    density=px_count/(length*rect.width)
    if density<density_thresh:     ##如果还达不到标准，则尝试减少seed的邻域半径
        return reduce_region_radius(rect,region,px_count,grad_mag,ll_angles,used,reg_angle,precision,p,density_thresh)
    ##达到这里，说明标准达到了
    return px_count,reg_angle,True

def rect_inter_low(cur_x,x1,y1,x2,y2):
    """
    给定当前x点，得到此列上与p1->p2线段的交点，
    取交点的y坐标作为下界
    """
    if x1==x2 and y1<y2:     return y1
    if x1==x2 and y1>y2:     return y2
    return y1+(cur_x-x1)*(y2-y1)/(x2-x1)

def rect_inter_hi(cur_x,x1,y1,x2,y2):
    """
    给定当前x点，得到此列上与p1->p2线段的交点，
    取交点的y坐标作为上界
    """
    if x1==x2 and y1<y2: return y2
    if x1==x2 and y1>y2: return y1
    return y1+(cur_x-x1)*(y2-y1)/(x2-x1)

def get_rect_pxs(Rect rect,double[:,:] ll_angles):
    """
    生成并返回给定rect对象包含的像素坐标列表
    """
    cdef:
        int px_cnt=0,alg_px_cnt=0
        int height,width,offset
        double dy,dx,ys,ye,cur_x,cur_y
        double tcorners[4][2]
        double corners[4][2]
    height,width=ll_angles.shape[:2]
    ##计算并调整矩形的角点
    dy=math.sin(rect.theta)
    dx=math.cos(rect.theta)
    tcorners[0][0]=rect.y1+dx*rect.width/2; tcorners[0][1]=rect.x1-dy*rect.width/2
    tcorners[1][0]=rect.y2+dx*rect.width/2; tcorners[1][1]=rect.x2-dy*rect.width/2
    tcorners[2][0]=rect.y2-dx*rect.width/2; tcorners[2][1]=rect.x2+dy*rect.width/2
    tcorners[3][0]=rect.y1-dx*rect.width/2; tcorners[3][1]=rect.x1+dy*rect.width/2
    if rect.x1<rect.x2 and rect.y1<=rect.y2 : offset=0
    elif rect.x1>=rect.x2 and rect.y1 <rect.y2: offset=1
    elif rect.x1>rect.x2 and rect.y1>=rect.y2: offset=2
    else: offset =3
    for i in range(4):
        corners[i][0]=tcorners[(offset+i)%4][0]
        corners[i][1]=tcorners[(offset+i)%4][1]
    cur_x=<int>math.ceil(corners[0][1])
    cur_y=<int>math.ceil(corners[0][0])
    ys=-1e+38; ye=-1e+38  #列的起始y坐标，终止y坐标
    while cur_x<corners[2][1] and cur_y>ye:
        ##计算得到当前x列的y坐标下界，上界
        if cur_x < corners[3][1]:
            ys=rect_inter_low(cur_x, corners[0][1],corners[0][0],corners[3][1],corners[3][0])
        else:
            ys=rect_inter_low(cur_x, corners[2][1],corners[2][0],corners[3][1],corners[3][0])
        if cur_x <corners[1][1]:
            ye=rect_inter_hi(cur_x,corners[0][1],corners[0][0],corners[1][1],corners[1][0])
        else:
            ye=rect_inter_hi(cur_x,corners[1][1],corners[1][0],corners[2][1],corners[2][0])
        ##添加在范围内的像素
        cur_y=math.ceil(ys)
        while cur_y<=ye:
            #是矩形中的点
            if cur_y >= 0 and cur_x >= 0 and cur_y < height and cur_x < width:
                px_cnt=px_cnt+1
                if isaligned(<int>cur_x,<int>cur_y, rect.theta,ll_angles,  rect.precision):
                    alg_px_cnt = alg_px_cnt + 1
            cur_y=cur_y+1
        cur_x=cur_x+1
    return px_cnt,alg_px_cnt

def log_gamma(double x):
    """
    计算对数的gamma函数
    """
    cdef double q[7]
    q[0]=75122.6331530
    q[1]=80916.6278952
    q[2]=36308.2951477
    q[3]=8687.24529705
    q[4]=1168.92649479
    q[5]=83.8676043424
    q[6]=2.50662827511
    if x>15.0:
        res=0.918938533204673 + (x-0.5)*math.log(x)-x \
            +0.5*x*math.log(x*math.sinh(1/x)+1/(810.0*x**6) )
    else:
        a=(x+0.5)*math.log(x+5.5)-(x+5.5)
        b=0
        for n in range(7):
            a=a-math.log(x+n)
            b=b+ q[n]*x**n
        res=a+math.log(b)
    return res

def rect_nfa(Rect rect,double[:,:] ll_angles,double logNT):
    """
    计算矩形的对数NFA值
    """
    #####计算矩形内总的有效像素个数，以及对齐点的个数
    cdef:
        int px_cnt,alg_px_cnt,n,k
        double log1Term,firstTerm,bin_tail,pterm,multiTerm,lastTerm
    px_cnt,alg_px_cnt=get_rect_pxs(rect,ll_angles)
    #####计算NFA
    n=px_cnt; k=alg_px_cnt
    if n==0 or k==0: return -logNT
    if n==k:
        return -logNT- n*math.log10(rect.p)
    log1Term=log_gamma(n+1)-log_gamma(k+1) \
              -log_gamma(n-k+1)+k*math.log(rect.p)+(n-k)*math.log(1-rect.p)
    firstTerm = math.exp(log1Term)
    ##在有些情况下计算第一项就足够了
    if abs(firstTerm)<2.2e-308:
        if k>n*rect.p:
            return -log1Term/math.log(10)-logNT
        else:
            return -logNT
    ##计算更多项
    bin_tail=firstTerm
    pterm=rect.p/(1-rect.p)
    multiTerm=0.0
    lastTerm=firstTerm
    for i in range(k+1,n+1):
        bin_term=(n-i+1)/i
        multiTerm=bin_term*pterm
        lastTerm=lastTerm*multiTerm
        bin_tail=bin_tail+lastTerm
        if bin_term<1.0:
            err=lastTerm*((1-multiTerm**(n-i+1))/(1-multiTerm)-1)
            if err<0.1*abs(-math.log10(bin_tail)-logNT)*bin_tail :
                break
    return -math.log10(bin_tail)-logNT


def rect_copy(Rect rect_src,Rect rect_dst):
    """
    将src中的信息，拷贝到dst中
    """
    rect_dst.x1=rect_src.x1;
    rect_dst.y1=rect_src.y1    #矩形主轴反方向上的端点
    rect_dst.x2=rect_src.x2
    rect_dst.y2=rect_src.y2    #矩形主轴正方向上的端点
    rect_dst.width=rect_src.width       #宽度
    rect_dst.cx= rect_src.cx     #矩形的中心点
    rect_dst.cy=rect_src.cy
    rect_dst.theta=rect_src.theta        #矩形的方向（弧度）
    rect_dst.precision=rect_src.precision
    rect_dst.p=rect_src.p
    rect_dst.dx=rect_src.dx
    rect_dst.dy=rect_src.dy


def rect_improve(Rect rect,Rect rect_bak,double[:,:] ll_angles,double logNT,double log_eps):
    """
    当矩形没有意义时，尝试矩形的一些变体来改进NFA值。
    返回值：细化的矩形，以及对数NFA值
    """
    cdef:
        int n
        double log_nfa,log_nfa_new
    log_nfa=rect_nfa(rect,ll_angles,logNT);
    if (log_nfa > log_eps):
        return log_nfa

    ##尝试细化精度
    rect_copy(rect,rect_bak)
    for n in range(5):
        rect_bak.p=rect_bak.p/2     #精度减小一半
        rect_bak.precision=rect_bak.p*PI
        log_nfa_new=rect_nfa(rect_bak,ll_angles,logNT)
        if log_nfa_new > log_nfa:
            log_nfa=log_nfa_new
            rect_copy(rect_bak,rect)
    if log_nfa > log_eps:
        return log_nfa

    ##尝试细化宽度
    rect_copy(rect,rect_bak)
    for n in range(5):
        if (rect_bak.width-0.5)>=0.5:
            rect_bak.width=rect_bak.width-0.5
            log_nfa_new=rect_nfa(rect_bak,ll_angles,logNT)
            if log_nfa_new > log_nfa:
                log_nfa=log_nfa_new
                rect_copy(rect_bak, rect)
    if log_nfa > log_eps:
        return log_nfa

    ##尝试细化一边的宽度
    rect_copy(rect,rect_bak)
    for n in range(5):
        if (rect_bak.width-0.5)>=0.5 :
            rect_bak.x1=rect_bak.x1 - rect_bak.dy*0.25
            rect_bak.y1=rect_bak.y1 + rect_bak.dx*0.25
            rect_bak.x2=rect_bak.x2 - rect_bak.dy*0.25
            rect_bak.y2=rect_bak.y2 + rect_bak.dx*0.25
            rect_bak.width = rect_bak.width-0.5
            log_nfa_new=rect_nfa(rect_bak,ll_angles,logNT)
            if log_nfa_new>log_nfa:
                log_nfa=log_nfa_new
                rect_copy(rect_bak, rect)
    if log_nfa >log_eps:
        return log_nfa
    ##尝试细化另一边的宽度
    rect_copy(rect,rect_bak)
    for n in range(5):
        if (rect_bak.width-0.5)>=0.5 :
            rect_bak.x1=rect_bak.x1 + rect_bak.dy*0.25
            rect_bak.y1=rect_bak.y1 - rect_bak.dx*0.25
            rect_bak.x2=rect_bak.x2 + rect_bak.dy*0.25
            rect_bak.y2=rect_bak.y2 - rect_bak.dx*0.25
            rect_bak.width = rect_bak.width-0.5
            log_nfa_new=rect_nfa(rect_bak,ll_angles,logNT)
            if log_nfa_new>log_nfa:
                log_nfa=log_nfa_new
                rect_copy(rect_bak, rect)
    if log_nfa >log_eps:
        return log_nfa

    ##尝试再次细化精度
    rect_copy(rect,rect_bak)
    for n in range(5):
        rect_bak.p=rect_bak.p/2     #精度减小一半
        rect_bak.precision=rect_bak.p*np.pi
        log_nfa_new=rect_nfa(rect_bak,ll_angles,logNT)
        if log_nfa_new > log_nfa:
            log_nfa=log_nfa_new
            rect_copy(rect_bak, rect)
    if log_nfa > log_eps:
        return log_nfa

    return log_nfa


def LineSegmentDetection(double[:,:] image, double scale, double sigma_scale,
                         double ang_thresh, double density_thresh,double quant,
                         int n_bins,double log_eps):
    """
    进行线段检测, 返加所有检测到的线段矩形对象列表
    """
    cdef:
        double precision,p,rho,logNT,reg_angle
        double[:,:] grad_mag
        double[:,:] ll_angles
        int[:,:] reg
        int[:,:] used
        int min_region_size,ls_count,i, px_count

    ##减采样
    scaled_image=gaussian_sampler(image,scale,sigma_scale)
    #scaled_image=self.gaussian_sampler(image,self.scale,self.sigma_scale)

    ##计算梯度幅度，角度，梯度伪序
    precision=PI*ang_thresh/180
    p=ang_thresh/180
    rho=quant/math.sin(precision)
    #grad_mag,ll_angles,list_p =self.levelline_angle(scaled_image,rho,self.n_bins);
    grad_mag, ll_angles, list_p = levelline_angle(scaled_image, rho, n_bins);
    height,width=ll_angles.shape[:2]
    used=np.zeros([height,width],dtype=np.int)
    ##logNT就是NFA公式中与图像宽高有关的项的对数
    logNT = 2.5 * ( math.log10( width ) + math.log10(height) ) +math.log10(11.0)
    min_region_size=<int>(-logNT/math.log10(p) )
    ##查找线段
    reg=np.zeros((height*width,2),dtype=np.int)
    ls_count=0
    rect_list=[]
    for i in range((height-1)*(width-1)):
        #px=tuple(px)
        x=list_p[i,1]; y=list_p[i,0]
        if used[y,x]==1 or ll_angles[y,x]==NOTDEF:
            continue
        #对当前像素执行区域生长算法
        #reg_angle,px_count=self.region_grow(px[1],px[0],ll_angles,reg,used,precision) #tested
        reg_angle,px_count=region_grow(x,y,ll_angles,reg,used, precision)
        if px_count <min_region_size:
            continue
        #为区域构建矩形
        rect=Rect()
        rect_bak=Rect()
        #rect=Rect()
        region2rect(rect,reg,px_count,reg_angle,precision,p,grad_mag)
        #self.region2rect(rect,reg,px_count,reg_angle,precision,p,grad_mag)
        #####refine
        #px_count,reg_angle,passed=self.refine(rect,reg,px_count,reg_angle,grad_mag,ll_angles,used,precision,p,self.density_thresh)
        px_count,reg_angle,passed=refine(rect,reg,px_count,reg_angle,grad_mag,ll_angles,used,precision,p,density_thresh)
        if not passed:
            continue
        #计算NFA值
        log_nfa=rect_improve(rect,rect_bak,ll_angles,logNT,log_eps)
        if log_nfa<=log_eps:
            continue
        ls_count=ls_count+1     ##找到并确认新的线段
        ##因为求梯度时是用的2x2掩膜，所以通过偏移校正
        rect.x1=rect.x1+0.5; rect.y1=rect.y1+0.5
        rect.x2=rect.x2+0.5; rect.y2=rect.y2+0.5
        if scale!=1.0:
            rect.x1=rect.x1/scale; rect.y1=rect.y1/scale
            rect.x2=rect.x2/scale; rect.y2=rect.y2/scale
            rect.width = rect.width/scale
        rect_list.append(rect)

    return rect_list